In [1]:
cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents

C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
In [2]:
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar

import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [3]:
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
In [24]:
from scipy.stats import boxcox
In [4]:
clima = gpd.read_file('climas_hemispheres')
In [5]:
clima.columns
Out[5]:
Index(['OBJECTID', 'GRIDCODE', 'is_north', 'SUM_2001_J', 'SUM_2001_F',
       'SUM_2001_M', 'SUM_2001_A', 'SUM_2001_1', 'SUM_2001_2', 'SUM_2001_3',
       ...
       'SUM_2017_3', 'SUM_2017_4', 'SUM_2017_S', 'SUM_2017_O', 'SUM_2017_N',
       'SUM_2017_D', 'SUM_Area', 'Shape_Leng', 'Shape_Area', 'geometry'],
      dtype='object', length=211)
In [7]:
area = clima['SUM_Area']*1e-6
In [9]:
data_columns = clima.columns[3:3+17*12]
In [10]:
data_columns
Out[10]:
Index(['SUM_2001_J', 'SUM_2001_F', 'SUM_2001_M', 'SUM_2001_A', 'SUM_2001_1',
       'SUM_2001_2', 'SUM_2001_3', 'SUM_2001_4', 'SUM_2001_S', 'SUM_2001_O',
       ...
       'SUM_2017_M', 'SUM_2017_A', 'SUM_2017_1', 'SUM_2017_2', 'SUM_2017_3',
       'SUM_2017_4', 'SUM_2017_S', 'SUM_2017_O', 'SUM_2017_N', 'SUM_2017_D'],
      dtype='object', length=204)
In [11]:
inter_cols = {}
for yr in range(17):
    inter_cols[yr] = data_columns[12*yr:12*yr +12]
In [12]:
intra_cols = {}
for mo in range(12):
    intra_cols[mo] = [data_columns[mo + yr *12  ] for yr in range(17)]
In [13]:
for yr in range(17):
    inter = clima[  inter_cols[yr]].sum(axis=1)/area
    label = 'inter_' + str(yr +2001)
    clima= clima.assign(label=inter)
    clima = clima.rename(columns={'label': label})
In [14]:
inter_columns = clima.columns[-17:]
In [15]:
inter_columns
Out[15]:
Index(['inter_2001', 'inter_2002', 'inter_2003', 'inter_2004', 'inter_2005',
       'inter_2006', 'inter_2007', 'inter_2008', 'inter_2009', 'inter_2010',
       'inter_2011', 'inter_2012', 'inter_2013', 'inter_2014', 'inter_2015',
       'inter_2016', 'inter_2017'],
      dtype='object')
In [16]:
inter_max = clima[inter_columns].max().max()
In [17]:
for mo in range(12):
    intra = clima[ intra_cols[mo]].sum(axis=1)/(17*area)
    label = 'intra_{:02}'.format(mo)
    clima = clima.assign(label=intra)
    clima = clima.rename(columns={'label': label})
In [18]:
intra_columns = clima.columns[-12:]
In [19]:
intra_columns
Out[19]:
Index(['intra_00', 'intra_01', 'intra_02', 'intra_03', 'intra_04', 'intra_05',
       'intra_06', 'intra_07', 'intra_08', 'intra_09', 'intra_10', 'intra_11'],
      dtype='object')
In [20]:
intra_max = clima[intra_columns].max().max()
In [21]:
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100,cmap='OrRd'):
    fig, ax = plt.subplots(figsize=(12, 8))

    ax.set_aspect('equal')
    ax = df.plot(ax =ax, column=column, cmap=cmap,vmin=vmin, vmax=vmax)
    ax.set_axis_off()
    
    norm = cls.Normalize(vmin, vmax)
    cmmapable = cm.ScalarMappable(norm, cmap)
    
    cmmapable._A = []
    cbaxes = inset_axes(ax, width="80%", height="1%", loc=8)
    cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal") 

    cb.set_label(title, fontsize=20, family='Times New Roman')
    cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')

    plt.show()
    fig.savefig(name, dpi=dpi, bbox_inches='tight')
In [22]:
# title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by climate in ' +mo_name + ' (2001-2017)'

inter-anual

In [25]:
lambdas = []
for yr in range(17):
    column= 'inter_{}'.format(2001 + yr)
    test_column, lbd = boxcox(clima[column] + 1e-25)
    lambdas.append(lbd)
lmbda=  sum(lambdas)/len(lambdas)
In [26]:
lambdas
Out[26]:
[0.14038468517156627,
 0.1444606660983037,
 0.15277534090253894,
 0.13317804176604167,
 0.14530905568227345,
 0.14102223123748542,
 0.14116689365003393,
 0.13685466780709593,
 0.14903029361538356,
 0.13685801384104096,
 0.13035492378545857,
 0.13746321513239707,
 0.14127730079045367,
 0.14911015883764694,
 0.14133583609742562,
 0.13494399135711255,
 0.14117937295672492]
In [27]:
maxs = []
mins = []
for yr in range(17):
    column= 'inter_{}'.format(2001 + yr)
    test_column = boxcox(clima[column] + 1e-25,lmbda=lmbda)
    maxs.append(test_column.max())
    mins.append(test_column.min())
In [28]:
box_min = min(mins)
box_delta =max(maxs) - box_min
In [29]:
for yr in range(17):
    column = 'inter_{}'.format(yr +2001)
    test_column = (boxcox(clima[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
    clima = clima.assign(test_column =test_column)
    yr_name = str(yr +2001)
    title = r'Burned area index (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +yr_name
    name = 'inter_climas_ns_index_{}'.format(yr +2001)
    draw_map(clima,'test_column',title=title,vmin=0,vmax=1,name=name)

C:\ProgramData\Anaconda3\envs\geopandas\lib\site-packages\matplotlib\font_manager.py:1339: UserWarning: findfont: Font family ['cursive'] not found. Falling back to DejaVu Sans
  (prop.get_family(), self.defaultFamily[fontext]))

intra-anual

In [30]:
lambdas = []
for mo in range(12):
    column= 'intra_{:02}'.format(mo)
    test_column, lbd = boxcox(clima[column] + 1e-25)
    lambdas.append(lbd)
lmbda=  sum(lambdas)/len(lambdas)
In [31]:
lambdas
Out[31]:
[0.09017671412617928,
 0.10625140833782082,
 0.11573628493147596,
 0.126511177066325,
 0.13046723613099198,
 0.1283594542227759,
 0.10990719723332583,
 0.1176504399986209,
 0.12326441044869325,
 0.11619317223433426,
 0.10079958089109725,
 0.09625699965587274]
In [32]:
maxs = []
mins = []
for mo in range(12):
    column= 'intra_{:02}'.format(mo)
    test_column = boxcox(clima[column] + 1e-25,lmbda=lmbda)
    maxs.append(test_column.max())
    mins.append(test_column.min())
In [33]:
box_min = min(mins)
box_delta =max(maxs) - box_min
In [34]:
for mo in range(12):
    column = 'intra_{:02}'.format(mo)
    test_column = (boxcox(clima[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
    clima = clima.assign(test_column =test_column)
    mo_name = calendar.month_name[mo + 1]
    title = r'Burned area index (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +mo_name + ' (2001-2017)'
    name = 'intra_climas_index_{:02}'.format(mo)
    draw_map(clima,'test_column',title=title,vmin=0,vmax=1,name=name)

total

In [35]:
total = clima[data_columns].sum(axis=1)/17

total = total/area
In [36]:
clima= clima.assign(total=total)
In [37]:
column = 'total'
test_column, lbd = boxcox(clima[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
clima = clima.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by climate in each hemisphere (2001-2017)' 
name = 'total_climas_nsparts_index'
draw_map(clima,'test_column',title=title,vmin=0,vmax=1,name=name, dpi=600)

Anomalias

Inter

In [38]:
inter_medias = clima[inter_columns].mean(axis=1)
In [39]:
clima[inter_columns].shape
Out[39]:
(49, 17)
In [40]:
anomalias = np.array(clima[inter_columns]) -  inter_medias.values.reshape(49,1)
In [44]:
bigger = abs(anomalias.max()) if abs(anomalias.max()) >= abs(anomalias.min()) else abs(anomalias.min())

vmax = bigger; vmin = -1 * bigger
In [45]:
for yr in range(17):
    test_column = anomalias[:,yr]
    clima = clima.assign(test_column =test_column)
    yr_name = str(yr +2001)
    
    title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +yr_name
    
    name = 'inter_anomalies_climas_ns_{}'.format(yr +2001)
    draw_map(clima,'test_column',title=title,name=name,cmap='seismic',vmin=vmin,vmax=vmax)

intra

In [47]:
intra_medias = clima[intra_columns].mean(axis=1)
In [48]:
anomalias = np.array(clima[intra_columns]) -  intra_medias.values.reshape(49,1)
In [49]:
bigger = abs(anomalias.max()) if abs(anomalias.max()) >= abs(anomalias.min()) else abs(anomalias.min())

vmax = bigger; vmin = -1 * bigger
In [50]:
for mo in range(12):
    test_column = anomalias[:,mo]
    clima= clima.assign(test_column =test_column)
    mo_name = calendar.month_name[mo + 1]
    
    title = r'Anomalies burned area pixels (#/$\mathrm{km}^{2}$) by climate in each hemisphere ' +mo_name + ' (2001-2017)'
    
    name = 'intra_anomalies_climate_ns_{:02}'.format(mo)
    draw_map(clima,'test_column',title=title,name=name,cmap='seismic',vmin=vmin,vmax=vmax, dpi =600)